rm(list = ls())


sink("log/log_replication_dangerous_liaisons.txt") # start log
begin.time<-Sys.time()



set.seed(123456789)

library(foreign)
library(bpNet)
library(ggplot2)
library(mcmcr)


## setwd('~/Desktop/PSRMreplicationLP/data/EmpData')
source('code/net_plot.R', chdir = TRUE)    ## load code to make figures

#Reading in trade matrices for each year, 
#converting them to matrices, and saving them as R-objects
#To do this, I first create strings for the filenames 
#to be read in (e.g., bin_trade1__1987.dta-bin_trade1_1999.dta) and 
#for the object names bin_trade1987-bin_trade1999) and then assign names to files.

N <- 126
TT <- 13

for(i in 1987:1999) {
filename = paste("data/EmpData/bin_trade1_", i,".dta", sep = "")	
  oname = paste("bin_trade_", i, sep = "")
  assign(oname, as.matrix(read.dta(filename)))
}


## create the W matrix 
W0 <- array(c(bin_trade_1987, bin_trade_1988,bin_trade_1989,
bin_trade_1990, bin_trade_1991, bin_trade_1992,bin_trade_1993,
bin_trade_1994,bin_trade_1995,bin_trade_1996,bin_trade_1997,
bin_trade_1998,bin_trade_1999),dim = c(126, 126, 13)) 

## set missing as 0 (and 10 as 0 because it appears on the diagnoal)
for (i in 1:13) {
  subw <- W0[,, i]
  subw[which(is.na(subw))] <- 0
  subw[which(subw == 10)] <- 0
  W0[,, i] <- subw 
}


## calculate indirect degree ratio 
indr <- function(w) {   ## w is a n*n binary matrix 
    
    N <- dim(w)[1]
    ## degree and indirect degree 
    d <- ind <- rep(NA, N)

    ## binary pos 
    pos1 <- sapply(1:N, function(i){return(which(w[i,] > 0))})
    pos2 <- sapply(1:N, function(i){return(which(w[i,] == 0))})

    ## calculate indirect degree 
    for (i in 1:N) {
        subpos <- pos2[[i]] ## have no direct link 
        ## remove itself 
        subpos <- setdiff(subpos, i)
        ## indirect link
        subind <- c()
        for (j in 1:length(subpos)) {
            subind <- c(subind, pos1[[subpos[j]]])
        }
        ## remove direct link 
        subind <- setdiff(subind, pos1[[i]])

        ## get measure 
        d[i] <- length(pos1[[i]])
        ind[i] <- length(subind)
    }

    return(ind / (ind + d))
}

## a matrix of indirect ratio N * TT
mind <- sapply(1:13, function(i){return(indr(W0[,,i]))})

## normalize the W matrix 
W <- W0
for (i in 1:N) {
    for (j in 1:TT) {
        subsum <- sum(W[i, , j])
        if (subsum > 0) {
            W[i, , j] <- W[i, , j] / subsum
        }
    }
}


#Creating the behavioral outcome matrix   
repress<-as.matrix(read.dta("data/EmpData/physint1.dta"))

# fill missings with unit-level average 
for (i in 1:N) {
  subm <- mean(repress[i, ], na.rm = 1)
  subna <- which(is.na(repress[i, ]))
  if (length(subna) > 0) {
    repress[i, subna] <- subm
  }
}


## convert to data frame 
id <- rep(1:N, each = TT)
time <- rep(1:TT, N)
Y <- c(t(repress))
indr <- c(t(mind))

data <- cbind.data.frame(id = id, time = time, Y = Y, indr = indr)

mcmc <- 35000
burnin <- 5000

out1.1 <- bpNet(data = data,   ## a data frame
                W = W,          ## an array N * N * TT, each slide is a w matrix for a given period
                index = c("id", "time"), ## id and time
                Yname = "Y",
                Xname = c("indr"),
                Zname = NULL,   ## covar has unit random effect
                Aname = NULL,   ##           time
                Contextual = NULL,  ## contectual covar
                Contextual.effect = "none",  ## contectual effect: "none", "unit", "time", "both"
                lagY = FALSE,                ## whether to incorporate lagged outcome
                force = "both",              ## two-way fe" "none", "unit", "time", "both"
                r = 10,
                flasso = 1,              ## factor selection
                rhoZ = as.matrix(rep(1, TT)),  ## time-varying covar that explains rho: T*k     
                rho_spec = "multilevel",                ## random rho: "multilevel", "rw", "ar1"
                niter = mcmc,
                burn = burnin)


## --------------------------------------- @@
## beta 
## simdata1
beta1 <- out1.1$beta

m.beta1 <- as.matrix(apply(beta1, 1, mean))
ci.beta1 <- t(apply(beta1, 1, quantile, c(0.025, 0.975)))

m1 <- cbind(m.beta1, ci.beta1)
m1 

cat("Model 1: Beta \n")
print(m1)
cat("\n")

## error variance 
pdf("plots/HR/hr_hist.pdf",width=18, height=9, pointsize = 20)
hist(as.mcmc(t(out1.1$sigma2)), 
     col="gray90", main = "", 
     xlab=expression(sigma^2))
dev.off()

## Figure A18: rho plot 
a <- 20
p1.1 <- rho_plot(out1.1, rep(0, TT), constant = FALSE)
p1.1 <- p1.1 +theme(axis.text.x = element_text( hjust = 1, size=a))
p1.1 <- p1.1 +theme(axis.text.y = element_text( hjust = 1, size=a))
ggsave('plots/HR/hr_rhot.pdf', p1.1, width = 18, height = 9)

## Figure A17: factors
p1 <- omegaPlot(out1.1, oxlim = c(-2, 2), 
                plotlength = 5, labelname = "Factor")
ggsave('plots/HR/hr_factor.pdf', p1, width = 18, height = 20)






out1.2 <- bpNet(data = data,   ## a data frame
                W = W,          ## an array N * N * TT, each slide is a w matrix for a given period
                index = c("id", "time"), ## id and time
                Yname = "Y",
                Xname = c("indr"),
                Zname = NULL,   ## covar has unit random effect
                Aname = NULL,   ##           time
                Contextual = NULL,  ## contectual covar
                Contextual.effect = "none",  ## contectual effect: "none", "unit", "time", "both"
                lagY = FALSE,                ## whether to incorporate lagged outcome
                force = "both",              ## two-way fe" "none", "unit", "time", "both"
                r = 10,
                flasso = 1,              ## factor selection
                rhoZ = as.matrix(rep(1, TT)),  ## time-varying covar that explains rho: T*k    
                rho_spec = "constant",                ## random rho: "multilevel", "rw", "ar1"
                niter = mcmc,
                burn = burnin)


beta2 <- out1.2$beta

m.beta2 <- as.matrix(apply(beta2, 1, mean))
ci.beta2 <- t(apply(beta2, 1, quantile, c(0.025, 0.975)))

m2 <- cbind(m.beta2, ci.beta2)
m2


cat("Model 2: Beta \n")
print(m2)
cat("\n")


## rho 
cat("Model 2: rho \n")
print(mean(out1.2$rho0))
cat("\n")



print(Sys.time()-begin.time) 
sink() # close log  



